%

tic;
sigma = 0.1;
n=3;

D = 1;

numtheta = 1000;
numr = 100;

epsilon = 2*pi*sigma/n;

N = 1:n;

theta0 = 0:2*pi/n:2*pi-2*pi/n;

r = 0:1/numr:1;
theta = 0:2*pi/numtheta:2*pi;

v = zeros(length(r),length(theta));
G = zeros(length(r),length(theta));
x = zeros(length(r),length(theta));
y = zeros(length(r),length(theta));

tstart = tic;
for i = 1:length(r)
    for j = 1:length(theta)
        for k = 1:n
            G(i,j) = G(i,j) - 1/pi * log(sqrt(1+r(i)^2-2*r(i)*cos(theta(j)-theta0(k)))) + r(i)^2/(4*pi) - 1/(8*pi);
        end
        v(i,j) = 1 / (D*n) * (-log(epsilon*n/2) + n/8 - pi*G(i,j));
        if v(i,j) < 0 
           v(i,j) = 0;
        end
    end
end

for i = 1:length(r)
    for j = 1:length(theta)
        x(i,j) = r(i)*cos(theta(j));
        y(i,j) = r(i)*sin(theta(j));
    end
end

tstop = toc
surf(x,y,v,'Edgecolor','none');
axis equal;